function [p,chi2stat] = coxgof(covariates,survivalTime,censorvec)


%%%%%
% This function computes the chi-square goodness-of-fit statistic for
% Cox regression models. This is an extension of the
% Hosmer-Lemeshow test to Cox regression.
% 
% The method is described in:
% R. B. D'Agostino and B. Nam. Evaluation of the performance of survival analysis models: descrimination and calibration measures. Handbook of Statistics, 23:1-25, 2003. 
%
% Inputs:
% covariates: 2D matrix containing covariates; rows are individuals and columns are covariates
% survivalTime - 1D vector containing survival times in days.
% censorvec - 1D boolean vector indicating right-censored data; 0 = not censored, 1 = right-censored. Should correspond to the elements in survivalTime.
%
% Written by Joon Lee, May. 12, 2011
%%%%%


[b,logl,H]=coxphfit(covariates,survivalTime,'censoring',censorvec,'baseline',0);
risk=exp(covariates*b);

[risk_sorted,idx]=sort(risk);
survivalTime_sorted=survivalTime(idx);
censorvec_sorted=censorvec(idx);

decileSize=round(length(risk)/10);
chi2stat=0;
O=zeros(10,1);
E=zeros(10,1);

% first 9 bins
for i=1:9
    first=(i-1)*decileSize+1;
    last=i*decileSize;
    [f,x]=ecdf(survivalTime_sorted(first:last),'censoring',censorvec_sorted(first:last));
    O(i)=f(end);
    E(i)=mean(1-exp(-H(end,2)).^risk_sorted(first:last));  
    chi2stat=chi2stat+(O(i)-E(i))^2/E(i)/(1-E(i));
end

% 10th bin (possibly with a different size than the other bins)
first=9*decileSize+1;
[f,x]=ecdf(survivalTime_sorted(first:end),'censoring',censorvec_sorted(first:end));
O(10)=f(end);
E(10)=mean(1-exp(-H(end,2)).^risk_sorted(first:end));
chi2stat=chi2stat+(O(10)-E(10))^2/E(10)/(1-E(10));

p=1-cdf('chi2',chi2stat,8);

